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We experimentally and numerically investigate the expansion of initially localized ultracold bosons 
in homogeneous one- and two-dimensional optical lattices. We find that both dimensionality and 
interaction strength crucially influence these non-equilibrium dynamics. While the atoms expand 
ballistically in all integrable limits, deviations from these limits dramatically suppress the expansion 
and lead to the appearance of almost bimodal cloud shapes, indicating diffusive dynamics in the 
center surrounded by ballistic wings. For strongly interacting bosons, we observe a dimensional 
crossover of the dynamics from ballistic in the one-dimensional hard-core case to diffusive in two 
dimensions, as well as a similar crossover when higher occupancies are introduced into the system. 



Non-equilibrium dynamics of strongly correlated 
many-body systems pose one of the most challenging 
problems for theoretical physics [1]. Especially in one di- 
mension, many fundamental questions concerning trans- 
port properties and relaxation dynamics in isolated sys- 
tems remain under active debate. These problems have 
attracted a renewed interest in recent years due to the 
advent of ultracold atomic gases. The ability to control 
various system parameters in real time has not only al- 
lowed quantum simulations of equilibrium properties of 
interacting many-body systems [2], but has also enabled 
experimental studies of quantum quenches [3-7] and par- 
ticle transport [8-12] in clean, well-controlled, and iso- 
lated systems. Here, we study the combined effects of in- 
teractions and dimensionality on the expansion dynamics 
of bosonic atoms in optical lattices. 

While interactions generally lead to diffusive trans- 
port in higher dimensions, the situation is more involved 
in one dimension (ID), where the phase space available 
for scattering can be severely limited. This was demon- 
strated, for example, by the experimental realization of 
a quantum Newton's cradle [5], showing that not all ID 
Bose gases thermalize (see also [13]). An intriguing phe- 
nomenon in ID is the existence of an exact mapping [14] 
from hard-core bosons on a lattice or a Tonks-Girardeau 
gas [15, 16] to non-interacting spinless fermions, demon- 
strating the integrability of these systems. Furthermore, 
this mapping establishes that the time evolution of the 
density distribution is identical for hard-core bosons and 
non-interacting fermions. As a consequence, hard-core 
bosons in ID expand ballistically and, asymptotically, 
undergo a dynamical fermionization during the expan- 
sion [17, 18]. In a transient regime, even initial ID 
Mott insulators with unity filling are predicted to be- 
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Figure 1. (Color) Experimental sequence and time evolu- 
tion during the expansion, (a) Sketch of the experimental 
sequence, (b)-(d) Experimental time evolution of line den- 
sity profiles during a ID expansion for various interaction 
strengths (each line is individually normalized), (e)-(g) Corre- 
sponding t-DMRG calculations for eight atoms, plotted using 
cubic interpolation. 



come coherent during the expansion and to dynamically 
form long-lived quasi-condensates at finite momenta [19- 
21]. In the presence of doubly occupied lattice sites 
(doublons) or even higher occupancies, the above map- 
ping is not applicable. The dynamics then become more 
involved and can include intriguing quantum distilla- 
tion effects, namely a demixing of doublons and single 
atoms [22, 23]. 

Several powerful theoretical methods have been used to 
study the expansion dynamics in ID, including the time- 
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dependent density matrix renormalization group method 
(t-DMRG) (see, e.g., [20, 22, 24]) and approaches based 
on the existence of exact solutions (see, e.g., [25-29]). 
For interacting two-dimensional (2D) systems, in con- 
trast, one needs to resort to approximate methods such 
as the time-dependent Gutz wilier ansatz, which predicts 
dynamical condensation even in 2D [30, 31]. 

In this work, we experimentally study the expansion 
of initially localized bosonic atoms in the lowest band 
of an optical lattice. We investigate how the expansion 
speed changes as a function of interaction strength and 
how it is affected by the dimensionality of the system. 
Furthermore, we identify the role of multiply occupied 
lattice sites in the system and compare our results to 
t-DMRG [32-34] calculations in the ID case. 

Experimental sequence - The experiment starts with 
a Bose-Einstein condensate of approximately 10^ bosonic 
^^K atoms in a three-beam optical dipole trap. The con- 
densate is loaded into a blue-detuned, three-dimensional 
optical lattice (lattice constant d = A/2, wavelength 
A = 736.7 nm) with a lattice depth of Vq = 33.0(5) E^. 
Here, E^ = / (2mA^) denotes the recoil energy, m the 
atomic mass, and h is Planck's constant. For suitable 
harmonic confinements, sufficiently strong repulsive in- 
teractions, and adiabatic loading, a large Mott insulating 
core with unity filling and a radius of (40 — 50) d is cre- 
ated in the center (see Fig. 1(a)). By employing a Fesh- 
bach resonance at a magnetic field of 402.50(3) G we can 
tune the interaction strength during loading and thereby 
control the amount of multiply occupied lattice sites. In 
the deep lattice, where tunneling is suppressed (tunnel- 
ing time Td = h/Jd ^ 58 ms, with the tunneling ampli- 
tude Jd and h = /i/(27r)), the atoms are held for a 20 ms 
dephasing period, during which any residual coherences 
between lattice sites are lost [35] and all atoms become 
localized to individual lattice sites. The resulting states 
after this loading procedure is a product of local Fock 

states, l^initiai) = ^ {biy^ |0) , Vi e {0, 1,2,...}, 

where P- is the creation operator for a boson on site i. 
This state is characterized by a fiat quasimomentum dis- 
tribution Uk = const, where k G [— Tr/d, Tr/d] denotes 
the quasimomentum. During the dephasing period, we 
change the magnetic field to set the desired interaction 
strength U/ J for the expansion. Due to the suppressed 
hopping during this part of the sequence, this field ramp 
does not alter the density distribution, i.e., the initial 
state prior to the expansion is identical for all interac- 
tions. The expansion is initiated by lowering the lat- 
tice depth along one or both horizontal directions (x, y) 
in 150 /is to a depth of 8.0 (1) E^ to induce tunneling 
with amplitudes Jx (r = h/Jx = 1.8 ms) and Jy between 
neighboring lattice sites along these directions. This is 
equivalent to a quantum quench from [// J oo to a 
finite U/J. Simultaneously, the strength of the dipole 
trap is reduced to a small but finite value that compen- 
sates the anti-confinement along the expansion direction 
created by the lattice beams (see [36] for supplementary 
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Figure 2. (Color) Core expansion velocity and dynamical gen- 
eration of higher occupancies, (a) Core expansion velocity 
Vc for experimental data in ID (black circles, lattice depth 
(8, 33, 33) along {x,y,z), Jy ~ 0) and 2D (blue circles, 
(8,8,33)^r, Jx = Jy) and t-DMRG calculations for iV = 10 
particles in ID (red triangles). Experimental error bars de- 
note the standard deviation of the linear fits. Inset: Vr calcu- 
lated by t-DMRG and extrapolated to infinite particle num- 
ber. Error bars are given by the uncertainty of the extrapo- 
lation [36]. (b) Higher occupancy, as measured by /h, versus 
expansion time in the experiment. For the points labeled 
"initial state", the measurement was performed directly after 
the dephasing period in the deep lattice [36]. (c) fh after an 
expansion time of t = 18 r. Error bars in (b) and (c) show 
standard deviations of averaging four data points. All lines 
are guides to the eye. 



information) . 

The dynamics in the resulting lattice can be described 
within the homogeneous Bose-Hubbard model: 

H = -JxY^ h% -JyY, blbj + ^ (n, - 1) . 



Here, U denotes the on-site interaction strength, hi = 
blbi^ and {i^j)x{y) indicates a summation over nearest 
neighbors along the x- (y-) direction. 

We monitor the in-situ density distribution of the ex- 
panding cloud using standard absorption imaging along 
the vertical axis. The recorded column densities are inte- 
grated over one direction and the resulting line densities 
are presented in Figs, l(b-d) as a function of the expan- 
sion time for the ID case. In both the non-interacting 
and the hard-core limits we expect a ballistic expansion 
which splits the cloud into a left- and right-moving por- 
tion [20, 37, 38], as can be seen in our numerical results 
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shown in Figs. l(e,g). While the sphtting can be clearly 
observed in the experimental data for the non-interacting 
case (Fig. 1(b)), the presence of a few multiply occupied 
lattice sites decreases its visibility in the strongly inter- 
acting case (Fig. 1(d)). 

Expansion velocities - To quantify the expansion 
dynamics we extract the half-width-at-half-maximum 
(HWHM) from the line density profiles and determine 
the core expansion velocities (Fig. 2(a)) via linear fits 
to the evolution of the HWHM at intermediate times [36] . 
In both ID and 2D, the maximum core expansion velocity 
occurs in the non-interacting limit, where the system ex- 
pands ballistically. Due to an exact dynamical symmetry 
of Hubbard models on bi-partite lattices, the expansion 
dynamics are independent of the sign of the interaction 
[39] and we therefore focus the discussion on the /7 > 
case. In 2D, increasing the interaction strength mono- 
tonically reduces the core expansion velocity until it es- 
sentially drops to zero. In ID, in contrast, a similar but 
much weaker suppression of the expansion velocity ex- 
tends only up to interaction strengths on the order of the 
bandwidth U ~ 47^^, while increases again for stronger 
interactions and eventually reaches values comparable to 
the non-interacting case. 

The same qualitative behavior is evident in the t- 
DMRG simulations for ten particles, shown as red 
triangles in Fig. 2(a). Since the numerically calcu- 
lated HWHM suffers from rather large finite-size ef- 
fects, we also present t-DMRG results for an al- 
ternative measure of the expansion velocity, namely 
= {d / dt) ^/ (t) — i^^(O), extracted from the radius 
R^{t) = {l/N)Y,.{ni{t)){i - i^fd^ (inset), where N 
is the particle number and iq denotes the central lat- 
tice site. It is more robust against finite-size effects 
and allows an extrapolation to infinite particle num- 
ber [36], and, in our setup, exhibits the same qualita- 
tive dependence on U. Moreover, at /7 = 0, has 
an intuitive physical interpretation, as it is in this case 
equal to the average expansion velocity The lat- 

ter is given by the initial quasimomentum distribution 
through = 1/ {Nh) ^^j,{dek/dk)'^nk , where Ck = 
—2Jcos{kd) denotes the tight-binding dispersion rela- 
tion. For the given initial state, where rik is flat, this 
results in Vg,^ = y/2{d/r)^ illustrated by the dashed line 
in the inset of Fig. 2(a). 

The fast expansion for strong interactions in ID is a 
consequence of the system entering into the hard-core bo- 
son regime, where, at = oo, it can be exactly mapped 
to non-interacting fermions, which expand ballistically 
with Vr = V^(<^/t) [36]. Even though hard-core bosons 
undergo collisions and their quasimomentum distribution 
changes over time [19, 20], the above mapping guarantees 
that the evolution of their density distribution is ballis- 
tic and identical to the non- interacting case. In other 
words, the conservation of the quasimomentum distribu- 
tion of the underlying non-interacting fermions severely 
constrains the scattering processes, thereby preventing 
the dynamics from becoming diffusive. 



Starting from the hard-core boson limit, the decrease 
of the expansion velocity towards smaller interactions can 
be qualitatively understood by considering the dynami- 
cal formation of doublons and higher occupancies. For 
U/J > 4, isolated doublons in ID can be thought of 
as heavy compound objects, propagating with typical ef- 
fective hopping matrix elements on the order of J'^/U 
[40] . While their formation is energetically suppressed at 
[// J 4, for smaller U the system can maximize its local 
entropy through the formation of doublons (and higher 
occupancies) during the early phase of the expansion (see 
Figs. 2(b,c)). Therefore, as U decreases, higher occupan- 
cies begin to form and the expansion velocity decreases. 
In addition, the possibility of creating higher occupan- 
cies increases the phase space available for scattering and 
therefore favors diffusive dynamics. For vanishing inter- 
actions, the scattering cross section approaches zero and 
the expansion becomes ballistic again with a large veloc- 
ity of Vr = V2{d/T). Therefore, there has to be a mini- 
mum of Vc at some intermediate /7, which turns out to be 
close to the critical U/J^ 3.3 for the ID Mott insulator 
to superfluid transition [41]. This is consistent with other 
studies of quantum quenches, which observe the fastest 
relaxation times close to the critical point [7, 42]. 

The buildup of higher occupancies during the initial 
expansion dynamics shown in Fig. 2(b) is monitored by 
comparing the number of atoms left after a parity pro- 
jection A^par with the total atom number A/totab yielding 
fh = (A^totai - A^par) /A^totai- In the abscucc of triply or 
higher occupied sites, /h measures the fraction of atoms 
on doubly occupied sites [36] . While the expansion starts 
from an initial state with essentially no higher occupancy, 
/h rises significantly over roughly the first half tunneling 
time. After this initial buildup, /h remains almost con- 
stant and changes only on the much slower timescale of 
the expansion (compare Fig. 2(c)). This initial fast relax- 
ation is purely local, as can be seen in t-DMRG calcula- 
tions comparing the relaxation timescale to the evolution 
of the system without opening the trap [36] . The forma- 
tion of higher occupancies is accompanied by changes in 
Uk and results in an increase of interaction energy and 
therefore a decrease in kinetic energy. The effect of the 
reduced kinetic energy (as measured by I'av) is, however, 
much smaller than the observed reduction of the expan- 
sion velocity [36]. We thus conclude that scattering pro- 
cesses during the expansion are mainly responsible for 
the slower expansion. 

1D-2D crossover - In Fig. 3 we analyze how the ex- 
pansion dynamics change when we gradually tune the 
dimensionality from a purely ID system towards a 2D 
geometry. This is implemented by varying the depth of 
the lattice along the ^/-direction and thereby the tun- 
neling ratio X — Jy/Jx foi" the expansion [43]. Upon 
increasing the expansion dynamics at strong interac- 
tions change fundamentally. Instead of the fast expansion 
observed in the ID case (Fig. 3(a)), the major fraction of 
the cloud simply remains in the center (Fig. 3(c)). More- 
over, the column density profiles shown in the insets of 
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Figure 3. (Color) 1D-2D Crossover, (a)-(c) Evolution of line 
density profiles for various tunneling ratios x — Jy/Jx and 
V jJx — 10. (d) Experimental core expansion velocity for 
various x- Lines are guides to the eye. Error bars denote 
the standard deviation of the linear fits. The insets show the 
column density at t ~ 36 r. 



Fig. 3(d) exhibit a characteristic bimodal structure. In 
the 2D case, this structure consists of a slowly expanding, 
round, diffusive core on top of a square-shaped ballistic 
background and can be seen for all moderate to strong 
interactions. In ID, on the other hand, a similar behavior 
is only visible for intermediate interaction strengths. 

In Fig. 3(d), we illustrate how the interaction depen- 
dence of changes as we go from a ID system with two 
integrable limits to a 2D system, where only the non- 
interacting case is integrable. The expansion speed in 
the non-interacting case is independent of x, since in this 
case the dynamics along the two lattice axes are separa- 
ble. For all values of x, the expansion speed initially 
decreases with increasing interactions. For small x^ the 
core expansion velocity increases again for strong inter- 
actions, whereas, for x > 0-5, it remains minimal. The 
behavior at large x, as well as the bimodal cloud shape, is 
analogous to the dynamics of strongly interacting lattice 
fermions in 2D, which were shown to be diffusive [11]. 
The square-shaped background consists of ballistically 
expanding atoms originating from the edge of the high 
density core, while collisions render the expansion diffu- 
sive for atoms inside the core. Such diffusive dynamics 
are consistent with the numerical observation that hard- 
core bosons in 2D thermalize [44]. Our experimental re- 
sults show a qualitative difference between the dynamics 
in ID and 2D in the strongly interacting regime, whereas 
theoretical studies using the time-dependent Gutzwiller 
ansatz predict a qualitatively similar behavior, indepen- 
dent of dimension [30, 31, 45]. Overall, we observe that, 
for interacting systems, the expansion along one direction 
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Figure 4. (Color) Effect of higher occupancies in the initial 
state. Core expansion velocity Vc in ID as a function of the 
higher occupancy, as measured by /h, in the initial state for 
various interaction strengths. Solid lines are guides to the eye. 
Vertical error bars show the standard deviation of the linear 
fits, horizontal error bars the standard deviation of averaging 
16 measurements. 



is suppressed by an increased tunneling along a trans- 
verse direction. This promotes the notion that increasing 
transverse tunneling enlarges the accessible phase-space 
for scattering processes and therefore favors diffusive dy- 
namics. 

Higher occupancies in the initial state - Figure 4 il- 
lustrates the effect of a random admixture of higher oc- 
cupancies in the initial state on the expansion dynam- 
ics. This admixture is created by loading the lattice at 
smaller interaction strength and higher densities, such 
that no clear Mott insulator will form. Nonetheless, the 
dephasing in the deep lattice remains effective, such that 
the initial state of the expansion can still be described 
as a product of local Fock states, but with higher oc- 
cupancies on some randomly chosen sites. While there 
is, as expected, no significant effect of multiply occu- 
pied sites in the non-interacting case, where each atom 
expands individually, already at U/Jx = 1 higher oc- 
cupancies in the initial state reduce the core expansion 
velocity. This reduction becomes most dramatic close to 
the hard-core limit (U / Jx = 10), where the originally 
high expansion velocity quickly approaches zero. In this 
limit, any higher occupancies are long-lived [46] and their 
small effective higher-order tunneling rate slows down the 
expansion [22, 24]. Furthermore, the presence of multi- 
ply occupied lattice sites in the strongly interacting limit 
can give rise to quantum distillation processes [22] and 
thereby the formation of a stable core of doubly occupied 
lattice sites [47]. 

Conclusion - Experimentally, we find the fastest ex- 
pansions near the exactly solvable limits of the Bose- 
Hubbard model, where additional conservation laws re- 
strict scattering such that diffusion is not possible. These 
are: (i) the non-interacting limit, irrespective of dimen- 
sion, and (ii) the case of infinitely strong interactions in 
ID, provided there are no higher occupancies in the ini- 
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tial state. Deviations from these cases, either by finite 
interactions, the crossover towards 2D, or an admixture 
of higher occupancies in the initial state, lead to a sub- 
stantial suppression of the expansion. In the case of the 
crossover to 2D at large /7/ J, the emergence of diffusive 
dynamics in the core is additionally signaled by the char- 
acteristic bimodal cloud shape previously observed in the 
fermionic case [11]. In ID at intermediate interactions or 
with initially multiply occupied lattice sites, both exper- 
imental and t-DMRG profiles suggest an almost bimodal 
structure here as well. Therefore, we conjecture that the 
common reason for the slow expansions seen in the ex- 
periments is the emergence of diffusive dynamics in the 



core region of the cloud. 
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Supplementary Material 



I. EXPERIMENTAL DETAILS 

A. Preparation of the initial state 

We prepare a condensate of approximately 10^ ^^K 
atoms in the |F = 1, ttif = 1) hyperfine state in a three- 
beam optical dipole trap with trap frequencies of ujx = 
27r X 52(2) Hz and ujy = 27r x 52(2) Hz along the hori- 
zontal {x,y) directions and cj^ = 27r x 119(8) Hz along 
the vertical {z) direction. For the experiments using an 
initial state of almost exclusively singly occupied sites, 
presented in Figs. 2 and 3 of the main text, the ini- 
tial scattering length is set to = 350(7) ao, where ao 
is Bohr's radius, by employing a Feshbach resonance at 
402.50(3) G [48]. We linearly ramp up the lattice poten- 
tial (lattice constant d = X/2 = 368.3 nm) in 8 ms to a 
depth of 20.0(3) Ej^. We freeze out the resulting density 
distribution by a second lattice ramp to 33.0(5) in 
1 ms. At the same time we turn off the vertical confine- 
ment by switching off the dipole trap beams along the 
horizontal axes. The atoms are now suspended against 
gravity solely by the deep vertical lattice, which remains 
at this depth during the rest of the experiment. The 
small tunneling amplitude of Jd = h x 2.7(2) Hz in the 
deep lattice, in combination with the effects of gravity, 
induces Bloch oscillations with an oscillation length of 
11(1) nm that is small compared to the lattice constant. 
The system is thereby effectively decoupled into inde- 
pendent 2D systems. The intensity of the vertical dipole 
trap beam is increased simultaneously with the loading of 
the lattice to compensate the increasing anti-confinement 
caused by the lattice beams. In the deep lattice, the com- 
bined trap frequencies along the horizontal directions are 
cja^ = 27r X 56(6) Hz and ujy = 27t x 50(5) Hz. The atoms 
are held in the deep 3D lattice for 20 ms while the mag- 
netic field is ramped to set the scattering length for the 
expansion. During this period, all coherences between 
lattice sites are lost and all atoms become localized to 
individual lattice sites [35]. The small tunneling rate 
ensures that atoms cannot redistribute during this hold 
time such that the resulting state is identical for all fi- 
nal interaction strengths. The expansion is initiated by 
ramping down the lattice depth along one or two horizon- 
tal directions and simultaneously adjusting the intensity 
of the remaining vertical dipole trap beam within 150 /is, 
a timescale that is slow enough to avoid excitations into 
higher bands but fast compared to the tunneling rate. 
Figure S5(a) presents a sketch of the intensity and field 
ramps employed in the experimental sequence. For the 
data with increased higher occupancy (Fig. 4), we re- 
duced the initial interaction strength and increased the 
initial trap frequencies, such that no clear Mott insulator 
could form. 
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Figure S5. Experimental sequence, optimization of homogene- 
ity and extraction of Vc . (a) Sketch of the experimental se- 
quence (not to scale) . (b) Width of Gaussian fit to line density 
profiles after an expansion time of t = 54t at U/Jx ~ 56 as 
a function of the trap frequency due to the vertical dipole 
trap beam, neglecting the anti-confinement due to the lattice 
beams. Insets show the column density distribution at the 
indicated points, (c) Extracted HWHM of expanding clouds 
in ID for various interaction strengths. The lines show the fit 
result for Vc and their extension indicates the fit range. 



B. Optimization of homogeneity 



The optical potentials during the expansion are pro- 
vided by the three blue-detuned optical lattices, with 
Gaussian waists of approximately 150 /im, and the red- 
detuned dipole trap beam along the vertical axis with 
approximately the same Gaussian waist. In the overlap 
region of these beams, the anti-confining potential due to 
the lattice beams can thus be compensated along the hor- 
izontal directions by the confining potential of the dipole 
trap beam. An exact compensation along both horizon- 
tal directions is possible for equal lattice depths along 
these directions. In all other cases, we optimize the com- 
pensation for the x-direction, along which we record the 
dynamics. To perform the optimization, we let the atoms 
expand in the lattice for a fixed expansion time and vari- 
ous intensities of the dipole trap beam and maximize the 
final size of the cloud (see Fig. S5(b)). Note that both 
confining as well as anti-confining potentials hinder the 
expansion [11]. 
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C. Determination of the core expansion velocity Vc 

To determine the core expansion velocity Vc along the 
x-direction, we first integrate the recorded in-situ column 
densities of the clouds along the transverse (y) direction 
to obtain line density profiles. For each of these pro- 
files, we determine the maximum density n-max- 

Start- 
ing from the approximate center of the profiles, as deter- 
mined by a Gaussian fit, we move outwards in both direc- 
tions and determine the positions where the density first 
drops to nmax/2, using linear interpolation between the 
points. Half of the distance between these two positions 
is recorded as the half-width-at-half-maximum (HWHM) 
of the cloud. Typical time evolutions of the HWHM are 
shown in Fig. S5(c). Even in the non-interacting cases, 
the HWHM does not significantly increase during the 
first few tunneling times of the evolution. The HWHM 
only grows significantly once the extension of the expand- 
ing single particle wavefunctions becomes comparable to 
the initial cloud size [28]. For very large clouds, the as- 
sumption of a homogeneous lattice with constant lattice 
depth is not valid anymore, because lattice and dipole 
trap beams have only a finite width, giving rise to resid- 
ual potentials. Thus, we apply a linear fit with offset 
to the HWHM evolution only in the time range from 
t ^ 14r to t ~ 42r. The slope of this fit is the core 
expansion velocity Vc. 

D. Error estimates for U and Jx 

Error bars for U/Jx and r are not given in Figs. 2 (a- 
c) and 3 (d) of the main text, as their statistical errors, 
caused by uncorrelated fluctuations, are much smaller 
than the width of the points. There are, however, sources 
of systematic errors for U and Jx. From measurements 
of the expansion speed of non-interacting bosons over 
several days, we estimate the long term fluctuations of 
our calibration of Jx to be on the order of 3%. The 
main contribution to the uncertainty in U originates from 
the uncertainty of the width of the Feshbach resonance 
that is used to calculate the scattering length at a 
given magnetic field. The resulting uncertainty of the 
scattering length ranges from 0.5 ao at a set value of ao 
to 3.8 ao at 188 ao- 



E. Detection of multiply occupied sites 

In order to determine the fraction of atoms on multi- 
ply occupied sites, we first freeze out the on-site num- 
ber distribution by ramping up the lattice in 50 fis to 
a depth of 33.0(5) Ej^ along all three axes. In the deep 
lattice, where tunneling is strongly suppressed (tunnel- 
ing time Td = h/J(i ^ 58 ms), we set the magnetic 
field within 10 ms to a fixed value of B ^ 400 G, where 
the scattering length is large. We then apply a near- 
resonant light pulse which is approximately 110 MHz 




Pulse duration fpuise(Ms) 

Figure S6. Detection of multiply occupied sites. Loss of atoms 
for off-resonant light pulses of varying duration. The fast loss 
(rh) originates from higher occupancies, the slow loss (ts) is 
caused by off-resonant excitations. 



red detuned relative to the high-field imaging transition 
from the |4^S'i/2,^/ = +3/2, mj = —1/2) state to the 
|4^P3/2,^/ = +3/2, mj = —3/2) state. On multiply oc- 
cupied sites, this near-resonant light pulse gives rise to 
a fast two-body loss process caused by light-assisted in- 
elastic scattering of atoms [49]. This loss acts as a par- 
ity projection of the on-site atom number and results 
in the loss of all atoms for even atom numbers and the 
loss of all but one atom for odd atom numbers. In 
Fig. S6, we present a typical decay curve of the total 
atom number in the presence of the near-resonant light 
for varying pulse durations. The initial fast loss (lifetime 
Th = 9.9(1.3) /is) stems from losses on multiply occupied 
sites. After the initial decay, we observe a much slower 
decay (ts = 399(13) jas) that is caused by off-resonant 
excitations of the remaining single atoms. We extract a 
measure of the higher occupancy by comparing the num- 
ber of atoms with (A/'puise) and without (A^totai) a near- 
resonant light pulse with a duration of tp^ise = 50 jas. In 
the presence of the near-resonant pulses, the parity pro- 
jection on multiply occupied sites has taken place and 
only atoms on singly occupied sites as well as the remain- 
ing atoms from sites with rji =3,5,... are left in the sys- 
tem. The measured atom number A^puise is then extrapo- 
lated to a pulse duration of /is using the measured slow 
decay time A^par = A^puise exp (tpuise/^s). We then cal- 
culate an approximate measure of the fraction of atoms 
on multiply occupied sites: /h = (A^totai - A^par) /A^totai- 
Note that, strictly speaking, /h is only a lower bound 
on the fraction of atoms on multiply occupied lattice 
sites, because A'par also contains one atom per site with 
r]i = 3,5,.... While being exact for singly and doubly 
occupied sites, the measured fractions will be systemati- 
cally too low whenever a significant amount of sites with 
> 3 is present in the system. However, significant 
amounts of sites with an occupation of > 3 contribute 
only in the weakly interacting regime (see the discussion 
of numerical results in Sec. II B 3). 
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II. TIME DEPENDENT DMRG SIMULATIONS 
FOR ID SYSTEMS 

We employ the adaptive time dependent density ma- 
trix renormalization group (t-DMRG) to carry out the 
time evolution [32-34, 50], using a Krylov-space based 
method to propagate the wavefunction in time [51, 52]. 
To ensure convergence of the time dependent wavefunc- 
tions we enforce a threshold on the discarded weight [53] 
per time-step St = 1/(16J) of 10~^. Typically, this cor- 
responds to using about 2000 states at the longest expan- 
sion times. Furthermore, we introduce a cutoff N]j = 3 in 
the number of bosons per site, unless stated otherwise. 
For the initial state considered in our numerical simu- 
lations, where the initial density is limited to (n^) = 2 
{hi = blbi), increasing the cutoff to A^b > 3 results in 
changes of only about 1% for most of the quantities con- 
sidered here. For /7/J < 2, however, the quasimomentum 
distributions are calculated with A/'b = since fluc- 
tuations in the local particle number are larger in this 
regime. In the two integrable limits, U = and = oo, 
we also use exact diagonalization to compute time de- 
pendent quantities. 

The time scales that can be accessed with t-DMRG 
are limited due to the growth of entanglement that is en- 
coded in the time-evolved wavefunction [50] . This growth 
depends on the type of non-equilibrium problem as well 
as on other factors such as particle number and system 
size. We want to reach timescales at which all particles 
participate in the expansion, which requires the total sys- 
tem size to be substantially larger (roughly by a factor of 
four in our simulations) than the extension of the region 
with a finite density in the initial state. Therefore, the 
particle number is restricted to N < 14. 



A. Initial state and measures for the expansion 
velocity 

The initial state for our simulations is a state with 
exactly 77^ = 0, 1 or 2 bosons per site: 

|V^initial)=n4^(^I)'^|0)- (SI) 

For the t-DMRG data shown in Figs. 1 and 2 of the main 
text and Fig. S7, we consider a state with = 1 on L 
adjacent sites and r]i = otherwise such that, in this 
case, the total particle number N = ^i{hi) equals L. 
For this particular initial state, the kinetic energy 

and the interaction energy 

^int = I ^{Mni - 1)> (S3) 



both vanish and the initial quasimomentum distribution 
nk = j^^e-^'^^'-"^\brb^) (S4) 

l,m 

is flat, Uk = const. The total energy E = £^kin + ^^int, 
often referred to as the release energy, is an integral of 
motion and its value is independent of U in our problem. 

In the experiment, there exists a harmonic trapping 
potential during the lattice loading, and hence also re- 
gions of lower average density, yet the particle numbers 
that we can run reliable simulations for are too small to 
account for this. 

We study several measures of the expansion velocity: 
First, the core expansion velocity Vc extracted from the 
time dependence of the HWHM, second, the expansion 
velocity = dR{t)/dt, where R{t) = ^R^{t) - R^{^), 
related to the time dependence of the radius of the cloud, 

RHt) = ^Y.^n,mi-i^?d\ (S5) 

i 

and third, we analyze the average velocity 




where = —2Jcos{kd) is the one-particle dispersion. 

We further compute the time dependent probability for 
multiply-occupied sites (i.e., sites with r]i > 1). Numer- 
ically, we compute the fraction of atoms on multiply 
occupied sites from 

= ]^ XlIZ ^ (^^'^) • (S 
i m=2 

where hm,i measures the probability of finding m bosons 
on site i, {r]i\hm,i\T]i) = Sr^i,m^ where It^^) are the local 
Fock states on site i. This should be compared to /h, the 
quantity that can be accessed in our experiment. Nu- 
merically, we compute it as 

/h = ^ Z(2("2,i) + 2(n3,i) + 4(n4,i) + ...)• (S8) 

i 

We have also studied initial states with a finite density 
of holes (i.e., some = in Eq. (SI)), surrounded by 
sites with rji = 1. We find that such single hole defects 
(in the absence of doubly occupied sites) do not influence 
the expansion velocity in the hard-core limit U = 00 and 
therefore we do not discuss these results here any further. 

B. Sudden expansion starting from initial states 
with exactly one boson per site 

1 . Time dependence of density profiles 

Figure S7 shows t-DMRG results for the time depen- 
dence of the density {hi{t)) (top row, (a-f)), the doublon 
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Figure S7. Time dependence of the density profile and of the probabilities for double and triple occupancy, (a)-(f) Density profile 
(ni(t)), (g)-(l) double occupancy (n2,i(t)), (m)-(r) triple occupancy {h3,i{t)), for the sudden expansion starting from a product 
of local Fock states with exactly one boson per site (A^ = 8, U / J — 0, 1, 2, 4, 20, oo). Panels (a), (g) and (m) show the results 
for free bosons, where the particles can form higher occupancies without energy cost. While the outer wavefronts in the time 
dependent particle density {ni{t)) are given by the maximum single-particle velocity of 2{d/T) for all [/, in panels (c) and (d) 
we observe a long-lived high density core significantly slowing down the expansion. The dynamics at U/ J — 20 (shown in (e), 
(k) and (q)) are very similar to the limit of hard-core bosons [//J = oo ((f), (1) and (r)). 



density (n2,i(t)) (middle row, (g-1)), and the density of 
triple occupancy {n^^i{t)) (bottom row, (m-r)). The den- 
sity profiles at [// J = 2 and 4 exhibit a bimodal struc- 
ture: fast, ballistic tails and a slowly expanding high- 
density core. 

The density of multiply occupied lattice sites is zero 
both in the initial state and at all times for /7 = oo. Af- 
ter opening the trap and for [// J < oo, multiply occupied 
sites are dynamically generated. A net production (re- 
sulting in an increase of i^h, Eq. (S7)) occurs due to initial 
relaxation dynamics following the quench to finite U / J . 
This has to be contrasted to expansions that start from 
a system that is in thermal equilibrium (compare, e.g., 
Refs. 11, 17, 19, and 37) for which a non-trivial time evo- 
lution is solely due to the quench of the trapping poten- 
tial to zero. Numerically, we observe that mostly double 
occupancies appear. Ai U = 0, multiply occupied sites 
do not have an effect on the expansion speed, whereas 
for U/J > 4 isolated doublons tunnel slower than sin- 
gle bosons since their effective hopping matrix element 
is ~ J^//7 at large U. In addition, we observe the ef- 
fect of quantum distillation [22] at U/J = 20, where the 
doublons move towards the center of the cloud and stay 
there up to the maximum simulation time. In the case of 
bosons (as compared to fermions), there is an attractive 
interaction between doublons, enhancing the stability of 
strings of doubly occupied sites over time [23, 47]. 



2. Cloud radius and expansion velocity 

We analyze the radius R{t) and the expansion velocity 
from data such as those shown in Fig. S7. Figure S8(a) 
shows the time dependence of the radius R{t) for = 10 



bosons and /7/J = 0,4,20,oo. In the two exactly solvable 
limits [7 = and oo, the radius increases linearly in time, 
i.e., R(t) = v^t with = 'yav(O), as expected for this bal- 
listic dynamics. In the U = oo case, this can be seen 
using the mapping to spinless fermions, whose quasimo- 
mentum distribution is time-independent n^(t) = const. 

At intermediate values of [/, we observe fast transient 
dynamics for t < r: up to roughly this point in time, the 
radius increases linearly with a slope that is independent 
of U. For t > r, it crosses over to a linear increase with a 
smaller slope, with the slope now depending on U (inset 
in Fig. SB (a)). The expansion velocity shown in Fig. 2 
of the main text and in Fig. S8(b) corresponds to this 
smaller slope. The deviation from a linear increase of 
the radius with a constant velocity is an indication for 
non-ballistic dynamics for small and intermediate values 
of U. 

In Fig. S8(b), we plot for finite particle numbers 
= 4 and = 10. We also extrapolated to A^ ^ oo 
by using data for A/" = 2, 4, 6, 8, 10 and a fitting function 
Vr{N) = Vr-\-a/N (inset). We only took into account data 
for the radius R{t) up to the same maximum time for all 
values of A^. We therefore cannot exclude a systematic 
error due to the limited timescales that can be reached 
in t-DMRG simulations (the larger A/", the shorter the 
accessible times). Qualitatively, we tend to overestimate 
Vr{N — ^ oo) in the vicinity of the minimum of v^. The 
result of this extrapolation is included in Fig. S8(b) as 
well and we conclude that the pronounced drop in the 
expansion velocity for < /7/ J < 10 is robust against 
finite-size effects. 

We stress that, in our case, the expansion velocity 
is not simply a measure of the release energy E = 
^kin + £^int7 which vauishes independently of U. This is in 
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Figure S8. Tzme dependence of the radius R of the expanding 
cloud, (a) R{t) for U/J = 0,4, 20,oo. Inset: Zoom-in on 
the transient dynamics for [//J = 0,4 and 8. (b) Expansion 
velocity Vr for = 4, 10 and extrapolated in 1/A^ to ^ oo. 
The dashed line indicates Vr = \^(d/r). Inset: Finite-size 
extrapolation for U / J — A and U/ J — 10. The solid line uses 
all data points, the dashed lines exclude either the smallest or 
largest N from the fits. The difference between the fit results 
of the solid line and the dashed lines determines the errors 
bars shown in the main graph. 




b Relaxation in trap 




Figure S9. Fraction of atoms on multiply occupied sites. We 
show both the exact fraction vi^ (solid lines) and the experi- 
mentally accessible measure /h (dashed lines), both computed 
with t-DMRG for the expansion from a state with rfi — 1 
(Eq. (SI)), for U/J = 0, 4, 10, 20, N = 10, and A^b = N. (a) 
Evolution during the expansion, (b) Evolution after a quench 
in [//J, without opening the trap. These data suggest that 
the transient time associated with the formation of higher oc- 
cupancies is given by t ~ 0.5 r such that this initial relaxation 
is purely local. 



contrast to previously studied expansions of interacting 
bosons [29, 54] or fermions [37] from an initial state that 
is in equilibrium. For instance, for the Tonks-Girardeau 
gas or weakly interacting Bose gases in a free space ex- 
pansion, in the asymptotic limit, cx \/E [29, 54]. 



3. Dynamical formation of multiply occupied lattice sites 

We argue that the fast transient dynamics evident in 
R{t) are due to local relaxation processes. It is instructive 
to consider two cases: (i) the sudden expansion under 
some value oi U / J realized in the experiment, and (ii) 
the time evolution without opening the trap, but after 
quenching to a finite value of U/ J < oo. (for theoretical 
studies on quantum quenches of the interaction strength 
in the Bose-Hubbard model, see, e.g., Refs. [55-63]). 

Figure S9(a) and (b) show the time dependence of the 
fraction of bosons on multiply occupied sites for U / J = 



0,4, 10,20 for cases (i) and (ii), respectively. We present 
both the total higher occupancy z^h (solid lines) and the 
experimentally accessible quantity /h (dashed lines) . For 
scenario (ii), the initial state Eq. (SI) with 77^ = 1 is 
not an eigenstate except ioi U = 00, and therefore the 
system explores phase space, resulting in a dynamical 
formation of multiply occupied lattice sites, i.e., z^h > 0. 
The fraction of atoms on multiply occupied sites is similar 
in cases (i) and (ii), corroborating the notion that the net 
production of higher occupancies is a local process, and 
therefore not a consequence of the expansion as such. 

Furthermore, the figure demonstrates that, for any 
U/J < 00, the system forms higher occupancies on a 
timescale of t ~ 0.5r. Both and /h saturate at [/- 
dependent values and, in the expanding case, slowly de- 
cay at larger times with small oscillations. Moreover, for 
U/J > 4, z^h ^ /h, indicating that in this regime only 
double occupancies are formed while higher occupancies 
are suppressed. Our t-DMRG results for /h are in qual- 
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Figure SIO. Time evolution of the quasimomentum distribu- 
tion nk(t). (a) Sudden expansion, (b) Relaxation without 
opening the trap. Both dynamics start from the state given 
in Eq. (SI) with ry^ = 1 at U/J = 1 for iV = 10 using Ni, = N. 



itative agreement with the experimental data presented 
in Fig. 2(b) of the main text concerning the timescales of 
the formation and the decay of higher occupancies. We 
ascribe quantitative differences to the presence of hole 
defects in the experiment and the inhomogeneity due to 
the harmonic trap in the experiment. 



4. Time dependence of the quasimomentum distribution 

The fast transient dynamics as opposed to the slower 
dynamics during the expansion can further be elucidated 
by considering the time dependence of the quasimomen- 
tum distribution n/c. In Fig. SlO(a) and (b), we display 
t-DMRG results for rikit) for the two cases (i) and (ii), 
respectively, with U / J = 1. In both scenarios, nk{t) de- 
velops a maximum at /c = on transient timescales t ^ r. 
This maximum remains stable in case (ii), while in case 
(i), the central peak slowly dissolves at later times during 
the expansion. 

From such data for nk{t) and for various values of /7, 
we calculate the average expansion velocity at time 
t = r as a function oiU / J and compare it to (Fig. Sll). 
Indeed, i'av(^ = t) of the expanding gas has a weak min- 
imum at U/J ~ 3 and increases for larger t//J, even- 
tually exceeding >/2((i/r). The latter, v^,^ > \/2((i/r), 
occurs because of the dynamical quasi-condensation [21] 
at large U where predominantly the quasimomenta at 
k = ±7r/{2d) become occupied [19, 20]. The minimum 
of ^av is present both for the actual expansion and the 
relaxation of the trapped gas without opening the trap. 
Therefore, we conclude that the decrease of the expansion 
velocity, measured through either Vc or v^, is partially 
due to the relaxation dynamics of the quasimomentum 
distribution at short times with a tendency of occupy- 
ing small momenta with a small velocity Vk <C 2{d/r) 
{vk = {l/h)dek/dk = (2d / t) sin {kd)) . However, as 
Fig. Sll clearly shows, <C '^av(^ = 't) at /7/J < 10, 
indicating that interactions during the expansion lead to 
an additional substantial drop of and 
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Figure Sll. Expansion velocity. Average expansion velocity 
a function of U/J for the expanding gas and without 
opening the trap, both at time t — r. We compare this to 
Vy. The dashed line indicates v — ^/2{d/r). All data are 
calculated for = 10 and v^v is calculated using A^b = A^ in 
both cases. 
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Figure S12. Sudden expansion in the presence of doublons in 
the initial state, (a), (b) Radius R{t) for U/ J = 10 and U / J = 
4. (c), (d) Time dependence of the HWHM for U/J = 10 and 
U / J — \. The data for /h 7^ were obtained by averaging 
over all distributions of doublons {r\i — 2) while keeping all 
other occupied sites in the initial state at ?7^ = 1. Dashed 
lines are linear fits to the last four data points that are used 
to obtain (N = 10 in all simulations). 



C. Expansion in the presence of doublons in the 
initial state 

In this section we investigate initial states that have a 
finite concentration of sites with an occupancy of 77^ = 2. 
Fixing the particle number to A" = 10, we generate all 
possible realizations with such defects for a given /h 7^ 
and calculate the averaged time dependent density. 

Figure S12 shows our results for the radius R{t) and 
the HWHM at U/J = 10 and U/J = 4, comparing the 
expansion from an initial state with rji = 1 only (/h = 0) 
to the ones with defects (/h 7^ 0). Already for the clean 
state there is a transient behavior in the HWHM before a 
linear increase in time sets in, which is consistent with the 
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experimental observations discussed in Sec. I C (compare 
Fig. S5(c)). A similar behavior of the HWHM was also 
seen in a theoretical study of the sudden expansion of a 
Tonks-Girardeau gas [28]. 

Upon adding doublons to the initial state, the slope 
of the radius R{t) decreases, as shown in Fig. S12(a) for 
U/J = 10. The effect is small since R{t) is dominated 
by the fast moving ballistic tails, which are unaffected 
by the presence of a few doublons (see Fig. S7). The 
HWHM is, however, much more sensitive to the presence 
of doublons in the initial state: already for /h = 0.2, its 



slope, is zero or slightly negative. We observe this 
dramatic dependence of Vc on /h for both U/ J = A and 
U/J = 10. These numerical results agree well with the 
experimental data for Vc = Vc{fh) shown in Fig. 4 of the 
main text. 

At smaller U/J ~ 1, our t-DMRG results do not show 
any strong effect of the presence of doublons in the initial 
state on either Vc or Vr, in contrast to the experimental 
results. We attribute this deviation to the different par- 
ticle numbers {N ~ 10 in t-DMRG simulations versus 
N ^ SO per tube in the experiment). 



